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ABSTRACT 

Context. 

Aims. Accurate measurement of gravitational shear from images of distant galaxies is one of the most direct ways of studying the distribution 
of mass in the universe. We describe a new implementation of a technique for measuring shear that is based on the shapelets formalism. 
Methods. The shapelets technique describes PSF and observed images in terms of Gauss-Hermite expansions (Gaussians times polynomials). 
It allows the various operations that a galaxy image undergoes before being registered in a camera (gravitational shear, PSF convolution, 
pixelation) to be modeled in a single formalism, so that intrinsic ellipticities can be derived in a single modeling step. 

Results. The resulting algorithm, and tests of it on idealized data as well as more realistic simulated images from the STEP project, are 
described. Results are very promising, with attained calibration accuracy better than four percent (1 percent for round PSFs) and PSF ellipticity 
correction better than a factor of 20. Residual calibration problems are discussed. 
Conclusions. 
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1. Introduction 

Weak gravitational lensing is recognized as a profitable way 
to study the dark matter distribution in the universe, and with 
a series of ever wider-field astronomical cameras coming on- 
line, very large weak lensing surveys are being planned and 
performed. The inevitable source of noise in weak lensing mea- 
surements is shape noise, caused by the diversity of projected 
galaxy shapes on the sky. To beat down the noise, very large 
numbers of galaxy ellipticities need to be averaged: thus new 
lensing surveys such as the CFHT Legacy Survey (Hoekstra et 
al.ED.05), the CTIO survey (Jarvis et al. l2003\ . the VIRMO S- 
Descart survey (van Waerbeke, Mellier & Hoekstra 2005} or 
the recently-approved Kilo-Degree Survey (KIDS) on the VLT 
Survey Telescope, will contain millions of background sources, 
which in principle should enable very accurate shear measure- 
ments. 

This averaging down of the shape noise through large num- 
ber statistics only makes sense if systematic errors can be con- 
trolled: the main one is still the blurring of the source images 
by atmosphere, telescope and detector pixels. The commonly- 
used Kaiser, Squires & Broadhurst ( 1995 ; henceforth KSB) and 
Luppino & Kaiser ( 1997 > methods provide recipes for correct- 
ing these effects, and are very successful. Nevertheless, they 
are based on an idealized model of the effect of point spread 
function (PSF) convolution on ellipticity, and it is possible to 
construct plausible PSFs that it fails to correct properly (e.g., 
Hoekstra et al. 1998 Appendix D). Therefore it seems unlikely 



that the KSB recipes will deliver the factor of 10 to 100 im- 
provement in fidelity that will be required to exploit the new 
surveys (Erben et al. 2001 1. 

A number of different approaches have been put forward to 
improve PSF correction (Kuijken 1 19991 Kaiser 127)001 Rhodes, 
Refregier & Groth[2UU0| Bernstein & Jarvis l^UOH Refregier & 
Bacon 2003 1 Mandelbaum et al. 2005 1. In this paper we present 
a new technique which combines elements from most of these. 

The paper is organized as follows: in section 2 shear and el- 
lipticity are defined, and the effect of one on the other. Section 
3 summarizes the shapelets formalism, and describes how a 
shapelet description of a source and its PSF can be used to gen- 
erate an ellipticity estimate that is useful for shear estimation. 
Section 4 substantiates the approach with idealized tests of the 
algorithm. In section 5 a software pipeline is presented that im- 
plements the full processing chain from an astronomical image 
to a shear estimate. Results of applying the pipeline to test data 
from the STEP project are shown in Section 6. In Section 7 we 
compare with other techniques, and Section 8 gives the conclu- 
sions. 

2. Preliminaries 

2.1. Shear and Distortion 

Following the usual practice, we parameterize the effect of 
weak gravitational lensing on a distant source in terms of 
a shear (71,72) and a convergence k: the distorted image 
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7i ense( j(x, y) is derived from the original I(x, y) via the transfor- 
mation Aensedfc y) = I(x',y'), where 
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(1) 



The first matrix in this equation (the distortion matrix) rep- 
resents the transformation from observed (x, y) to undistorted 
(x',y') coordinates. 

Without knowledge of the intrinsic source size, only the 
distortion (81,82) = (7i>72)/(l - «), which affects the shape of 
the source, can be measured (Schneider & Seitz ll995l . 

2.2. Ellipticity 

We define the ellipticity of an object's image / as follows: 

Let 7 e ii be the model image with constant-ellipticity 
isophotes that best approximates /. Then the ellipticity (e\,e2) 
of / is defined such that a distortion of (— e\, -e 2 ) makes / e u 
circular. 

The major axis position angle (p and the axis ratio q of an 
elliptical source are simply related to (e\, e 2 ): 



1 = 



1 -, 

IT 



and 



e\ — e cos . 



e2 — e sin 2<p (2) 



This definition is similar to the one adopted by Bernstein & 



Jarvis (2002 henceforth BJ02), but does not explicitly force a 
fit to an elliptical Gaussian. 

As discussed by BJ02, expressing the shapes of objects in 
terms of distortions e has a practical advantage: in this formu- 
lation it is simple to calculate the response of object shapes to 
small distortions. An elliptical source with ellipticity (e\,ei) 
that is sheared by a small amount (g\,g2) can be viewed as a 
circular source that is sheared twice, first by e, and then by gi, 
giving a combined distortion matrix 

/ 1 -e\ -«2 Ul-gi -82 
\ -<? 2 1 +e\ J\ -g 2 1 + gi 

_ / 1 - e\ ~ 81 + e \g\ + e 2g2 -e 2 - g2 + e\g2 - e 2 gi 



-e 2 - g2 ~ e\g2 + e 2 gi 1 + ei + g\ + eigi + e 2 gi 



(3) 



This matrix is no longer a pure distortion matrix (which would 
have to be symmetric and of trace 2), but some algebra shows 
that this matrix can be decomposed into a magnification, a ro- 
tation and a distortion: 

/ 1 -<?i -e 2 \( 1-gi -g 2 
\ -e 2 1 + e\ J\ -g 2 1 + gi 

1 R \ I 1 — e\ — 8\ —e 2 — 5 2 
—e 2 - S 2 1 + e\ + 61 



-R 1 

where, to first order in the distortion gi, 



(4) 
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e\g\ + e 2 g2 
e\g2 ~ e 2 g\ 
(l-el-ehg! 



5 2 = (1 -e\-e l 2 )g 2 



(5) 
(6) 
(7) 
(8) 



Thus the action of a small distortion g on an elliptical source 
with ellipticity e, (according to the definition in eq.0 is equiva- 
lent to acting on a circular source with, successively, a magnifi- 
cation, a rotation (neither of which affect the shape of a circular 
source), and a distortion by e, + 5,. 

Assuming now that we have an ensemble of elliptical 
sources, of random orientations, so that before distortion (e, ) = 
and {e 2 .) + {el) = (e 2 ), the average ellipticity of the population 



after a distortion (g\. 



g 2 ) is simply 



(9) 



3. Shapelets 

The shapelets basis is described in Refregier ( 2003 ). It consists 
of the two-dimensional Cartesian Gauss-Hermite functions, fa- 
mous as the energy eigenstates of the 2-D quantum harmonic 
oscillator: 

B ah (x,y) = r i /3- 1 e- [(l ^ )2+(v -- ,v)2l/2/s2 // fl (x/ y 6)//*(y/^). (10) 

Here (x, y) are coordinates on the image plane, x c and y c are the 
center of the expansion, B ah is the basis function of order (a, b), 
H a is the Hermite polynomial of order a and /3 is a scale radius. 

k" is a normalization constant chosen so that (B ab ) integrates 
to one. 

Shapelets are a convenient basis set for describing astro- 
nomical images because of the compact way in which vari- 
ous operators (translation, magnification, rotation, shear) can 
be expressed as matrices that act on the shapelet coefficients. 
Shapelets have a free scale radius B (the size of the Gaussian 
core of the functions), and R03 shows how the coefficients 
transform under change of B, and how to convolve objects with 
different scale radii. 

To avoid introducing a preferred direction, the expansion 
should be truncated in combined order N - a + b, not in a 
or b separately. (A basis set truncated in is complete under 
rotation. Effectively, such a truncation describes an image as a 
product of a circular Gaussian with an inhomogeneous polyno- 
mial in x and y of order Af. Rotation of such an image will mix 
the x'y-i terms at constant i + j.) 

The reason for choosing shapelets as a formalism for weak 
lensing analysis is its ability to describe the main operations 
that a galaxy image undergoes before it is registered at a tele- 
scope focal plane: in reverse order, pixelation, convolution with 
a PSF, and distortion. 

In this paper we concentrate on well-sampled (PSF FWHM 
at least 3^1 pixels), ground-based seeing-limited images. It re- 
mains to be seen to what extent diffraction-limited PSFs, and 
undersampling, can be handled with this formalism. 

3.1. Pixelation 

An image that is registered on a CCD is pixelated: the flux 
on the surface of the detector is read out in binned form. 
Mathematically, the flux is first boxcar smoothed (i.e., con- 
volved with a pixel), and then sampled at a spatial frequency 
of once per pixel. Therefore, if we fit a shapelet expansion to 
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the binned image I(k, I) as a linear superposition of the basis 
functions B" b : 

N N-a 

I(kJ) = Y,J] s <>bB ab (k- Xc ,l-y c ) (11) 

a=0 6=0 

then the shapelet coefficients s ab describe the pixel-convolved 
image directly. If a similar fit is made to a PSF image that is pix- 
elated in the same way, then the intrinsic, pre-seeing image is 
exactly the deconvolution of the two shapelet expansions (apart 
from the effects of noise, undersampling, and truncation of the 
shapelet expansions). 

3.2. PSF convolution 

Convolution with the PSF can be expressed as multiplication 
with a PSF matrix P fl] 2, lfl2 i, 2 O8in,y0out): if the shapelet coefficients 
of the PSF are p ab , and those of a model source are m ab , then 
their convolution is 

(la,bPabBfJ®(Z a , b m ab B£) 

Q3in,/3 udm aib2 )B2 b \ (12) 

The subscripts on B ab identify the scale radius B, which can 
be different for the three shapelet series involved: those for the 
PSF, input model source and output result of the convolution. 
P(Bi n ,B out ) convolves a shapelet expansion with scale radius B m 
with the PSF, resulting in a shapelet expansion with scale radius 
B oul . The coefficients of the PSF matrix are 

1 a\b\a 1 b 1 WmiPaxA) — / ; ^a^aj y ~-b l b 2 b } V<*s<*i V 1 -^ 
ay, by 

Here the elemental convolution matrix (f 3 ^ 1 expresses the 
convolution of two one-dimensional shapelets of scales B\ and 
B 2 as a new shapelet series with scale B3. A recurrence relation 
for C„ m i is given in R03. 

3.3. Ellipticity from shapelets 

Given a PSF and a PSF-convolved source, both expressed as 
shapelet series, we determine the ellipticity of the source by 
modeling it as a PSF-convolved, distorted circular source of 
arbitrary radial brightness profile. This approach is similar to 
the one described in Kuijken ( 1999), but it is more effective 
when expressed in terms of shapelets. 

A circular source of arbitrary radial brightness profile can 
be written as a series of circular shapelets C" of the form cqC° + 
c 2 C 2 + C4C 4 + . . ., where the c„ are free coefficients. The C" (see 
Appendix) are normalized to have unit integral over x, y, so c n 
gives the total counts in each component. After distortion of 
such a source, it becomes an elliptical source which, to leading 
order in ellipticity <?, can be written as 

(1 + eiSi + e 2 S 2 )(c Q C° + c 2 C 2 + c 4 C 4 + . . .) (14) 

where Si are the first-order shear operators (see R03). 

This machinery allows us to write the model for the ob- 
served source as 

P ■ (1 + eiSi + e 2 S 2 )(c C° + c 2 C 2 + c 4 C 4 + ...), (15) 



expressed as a set of shapelet coefficients that depend on the e, 
and the c„. Fitting this model to the shapelet coefficients of the 
observed source (with their errors) yields the best-fit ellipticity 
(e\, e 2 ) and the associated errors. 

To improve the accuracy, we make two modifications: we 
add centroid error parameters, and we only fit the model to a 
subset of the shapelet expansions. The centroid error parame- 
ters are included to allow for a mismatch between the center 
of the object and its shapelet expansion. If the PSF and/or the 
galaxy have some lopsidedness to them, the center of their best- 
fit shapelet expansion may not be at the flux-weighted center 
of the source (since our centering technique simply requires 
the 01 and 10 components to be exactly zero). Hence the cen- 
troid may move under convolution, which would spoil the fit. 
To guard against this, instead of fitting Eq.^]we fit a model of 
the form 

P-(l+eiSi+e 2 S2+diTi+d2T2)(coC +C2C 2 +. . .+c Nc C Nc ) (16) 

instead. The free translation terms djT\ in the model are ex- 
pressed in terms of the shapelet operators Ti. 

The second modification is made to contain truncation er- 
rors. While the shapelet basis is complete, and hence can de- 
scribe any source given enough terms, in practice the fact that 
the source is only sampled in a finite set of pixels means that 
the expansion needs to be truncated. Hence, except in very spe- 
cial cases, the PSF and galaxy are not described perfectly by a 
truncated shapelet series. The missing information propagates 
through the analysis, and is a source of systematic error in the 
PSF convolution (as some PSF terms may be missing) and in 
the calculation of the action of shear (since shearing high-order 
shapelets generates also lower-order terms). 

Truncation effects can be seen most clearly if we re-express 
the shapelets in polar (r, 9) coordinates (they become Gaussians 
times Laguerre polynomials of r — see BJ02, R03, Massey & 
Refregier 2005 1. Polar shapelets are combinations of Cartesian 
shapelets B ab of the same order N = a + b, whose angular 
dependence is a pure sine or cosine of md, for angular order m. 
The order of the Laguerre polynomial is N, with N > m and 
N + m even. 

The translation and shear operators, when applied to a polar 
shapelet of order (N, m), generate terms at order (N ± 1 , m + 1) 
and (N ±2,m + 2), respectively. If we truncate the shapelet 
series of the best-fit circular model for the pre-seeing, pre-shear 
galaxy at order N c , then to be consistent the m = 1 and m = 2 
series must be truncated at order N c — 1 and N c -2, respectively. 

Note that we are never completely safe from truncation 
effects: in particular, complex PSFs can in principle mix co- 
efficients of all orders. The problem is minimized, but not 
completely eliminated, by adopting suitable scale radii so that 
the amount of information carried in the high-order coeffi- 
cients is small. We further include a 'safety margin' by setting 
N c - N — 2, in case the highest-order shapelet coefficients are 
affected by PSF structure at even higher order. The scheme is 
illustrated for the typical case of = 8 in Fig.^ The highest- 
order polar shapelet coefficients that should be included in the 
fit are (N c , 0), (JV C - 1, ±1) and (JV e - 2, +2). 

In summary, from the shapelet series for each source, up to 
Cartesian order A^ max , we form a truncated polar shapelet series 



4 



Konrad Kuijken: Shears from shapelets 




8.8 
8,6 
8,4 
8,2 
8.0 
8,-2 
8,-4 



8,-6 



8,-8 



Cartesian (ab) 



Polar (Nm) 



Fig. 1. The information used in the ellipticity determination, for shapelet expansion to order N = 8. On the left, the Cartesian 
shapelet coefficients that are fitted to describe a source and its associated PSF. On the right, the same information has been 
rearranged into a polar shapelet expansion (the two may be transformed into one another by appropriate mixing of the terms at 
order N = a + b). The heavy line indicates the polar shapelets from which the ellipticity is estimated by a fit to eq.[H)] Note the 
safety margin at order N — 1 and N, and the fact that only azimuthal orders between -2 and +2 are fitted. 



including terms of order (m = 0, N = 0, 2, . . . , A^ max - 2), (m = 
l,N= l,3,...,iV m ax-3)and (m = 2,N = 2,4, . . ..AW - 4). 
The effect of the shear, translation and PSF operators on circu- 
lar basis functions up to order Af max - 2 is then calculated up to 
order A^x, and the result likewise converted to polar shapelets 
up to order N max - 2. 

Least-squares fitting the model to each source yields an el- 
lipticity estimate (e\,ei), expressed as the shear that needs to 
be applied to a circular source to fit the object optimally. 

Performing the least-squares fit is straightforward to do nu- 
merically. x 1 is a fourth-order polynomial in the fit parame- 
ters {co, C2, . . . , cn-2, e\,ei,d\, a^}, and the minimum can typi- 
cally be found in a few Levenberg-Marquardt iterations (Press 
et al. [T986V 

The errors on the shapelet coefficients for each source can 
be derived from the photon noise, and these can be propagated 
through in the^- 2 function. The second partial derivatives of x 2 
at the best fit give the inverse covariance matrix, which can be 
inverted to show the variance and covariances between the fit 
parameters. This results in proper error estimates on all param- 
eters, in particular on e\ and e2- In practice the errors on e, are 
only weakly correlated with those on di and c„. 



4. Tests 

We now describe some elementary tests of this approach. 



4.1. Test of ellipticity estimates 

For small ellipticity e the linear shear operator provides a good 
description of the action of the shear, but for larger e higher- 
order corrections come into play, and these corrections de- 
pend on the radial brightness profile. We have calibrated these 
corrections empirically using a set of model sources that fol- 



low a Sersic (1968 1 distribution of brightness, of index 1-4. 
Each source was sheared by varying amounts and encoded into 
shapelets using a range of different scale radii. The first-order 
ellipticity estimate ei st derived by fitting a model of ea.H4lwas 
then compared to the true ellipticity (see Fig. |2}. For small e 
the correct ellipticity is recovered, but at larger e the discrep- 
ancy grows. Fitting the residuals versus the radial profile shape 
parameters cq, Ci, c\, it turns out that the true ellipticity <? tlue 
can be derived from the fitted lst-order estimate by applying 
the correction factor 



-true 

eist 



1 



,2 
- 1st 



-0.41 + 



0.085c 2 + 0.63c 4 



(17) 



The formula is valid to better than 1% accuracy for e < 0.7, 
corresponding to axis ratios of nearly 6:1. 

Below we will, in fact, NOT apply this correction, because 
the errors on c , C2 and C4 are typically so high, and (in the case 
of small, barely resolved sources) correlated, that the correc- 
tion factor cannot be determined accurately. Fortunately most 
galaxies in the sky have ellipticity below 0.3, where the correc- 
tion is below 3%. If necessary, the accuracy could be further 
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c 2 /c o'/c 

Fig. 2. Empirical calibration of the post-linear correction to the 
measurement of e. Top left: fractional error on derived lst- 
order ellipticity e for a range of different scale parameters and 
Sersic indices. Top right: fractional error of the corrected el- 
lipticity. Bottom left: the coverage of the C2/C0, C4/C0 plane by 
the models. The four groups of points correspond, from top to 
bottom, to Sersic index 4, 3, 2 and 1. The horizontal spread is 
mostly a consequence of using different scale radii (3. Bottom 
Right: residuals on lst-order ellipticity vs. the linear combina- 
tion of c' Icq = (0.085C2 + 0.63c4)/co for input ellipticity 0.4, 
showing that this parameter drives the scatter. 




-40 -20 20 

X 



Fig. 3. Different shapelet fits to a Moffat PSF (1 + ar 2 y 2 with 
FWHM of 12 pixels. Top panel: the solid line shows a cross- 
section of an 8th-order fit, using the dispersion of the best-fit 
Gaussian as scale radius. The dotted line shows the correspond- 
ing fit for a scale radius that is a factor of 1.3 times larger. The 
lower panel shows the PSF multiplied by radius 3 , in order to ac- 
centuate the residuals in the wings. Dots are the actual (monte- 
carlo realized) pixelated PSFs used in the simulations of |4] 
Note how the increased scale radius makes for a much better fit 
in the outer regions, without a serious degradation in the core 
of the PSF. 



represents a good compromise, though its exact value is not 
critical. 



increased by evaluating the effect of shear on the shapelet basis 
functions to higher order in eq.1161 

4.2. Choice of Scale Radius 

A truncated shapelet expansion can only describe deviations 
from a Gaussian over a particular range of spatial scales (which 
widens with order N). For ellipticity determinations the outer 
parts of galaxy and PSF images are most important (e.g., the 
classic second moments depend on the 3rd moment of the ra- 
dial profile), so there is some advantage to taking as large a 
scale radius as possible. On the other hand, this radius should 
not be so large that the inner structure of the source cannot be 
resolved. 

We have found that, for shapelet expansions up to order 
N = 8 or higher, taking a scale radius which is 1.3 times the 
dispersion of the best-fit Gaussian works well for a range of 
model PSFs. Fig.[3]shows an example for a Moffat PSF with 
index 2: both the core and the wings can be fitted adequately 
with this choice of 

The impact of the choice of scale radius on ellipticity mea- 
surements is shown in Fig.|4j for typical images. The factor 1.3 



4.3. Test of PSF correction 

Correcting for the effects of PSF convolution (dilution of ellip- 
ticity by a round PSF, or biasing of ellipticity by an elongated 
one) is the most critical part of weak lensing. The following 
tests show how well the shapelets technique can do this. 

For the tests we generate simulated, high signal-to-noise 
(S/N) sources. To be sure that pixelation effects are taken into 
account properly, and that the PSF convolutions are done ac- 
curately, a brute-force technique is used: each source is built 
of individual 'photons' that are drawn from a 2-D Sersic dis- 
tribution, sheared if required, then have a 'PSF' displacement 
added to them, and finally are added to the pixel in which they 
fall. We use 10 million photons per source, which gives effec- 
tively noise-free images (S/N^;1000). 

We ran three sets of tests. In each case we explore galaxies 
with Sersic laws f(r) ~ e\p(-kr l ^" s ) with indices n s = 0.5 
(Gaussian), 1 (exponential), 2, 3 and 4 (de Vaucouleur), and 
use PSFs with a Moffat profile (1 + ar 2 )' 13 '" of index j3,„ = 2, 3 
and 9 (nearly Gaussian). All galaxies are scaled to an intrinsic 
effective radius of 4 pixels, and the PSF FWHM's range from 
4 to 12 pixels. We use shapelet expansions to order N - % and 
12, and scale radii of 1.3 times the dispersion of the best-fit 
Gaussian. 
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Fig. 5. Fractional error in recovering a 10% shear, using round PSFs. Each panel represents a different Moffat PSF; the rightmost 
panels are very nearly Gaussian. The simulated galaxies have effective radii of 4 pixels. Top row: 8th-order shapelets; bottom 
row: 12th-order shapelets. 



First, to test the 'shear calibration' factor, we check how 
well we can recover the shear of a galaxy that is sheared by 
g\ =0.1 and convolved with a round PSF. Fig. [5] shows the 
result: any calibration error is at the sub-percent level; the worst 
results are obtained for the most non-Gaussian PSFs (n s = 4). 
The noise in the curves suggests that we are also limited by the 
accuracy of the Monte-Carlo simulations of the galaxies, and 
by the numerics of the software implementation of the method. 

The second test shows to what extent PSF ellipticity pol- 
lutes the shear estimate. The input PSFs were given an elipticity 
e\ =0.1 (axis ratio 0.82), and convolved with round galaxies. 
The recovered e is shown in Fig. [6] The residual effect is at 
most half a percent (worst case), which represents a correction 
of the PSF ellipticity by a factor of better than 20. The best re- 
sults are obtained for lowish Sersic indices (below 2) and for 
PSFs with not too extended wings (Moffat index 3 or higher). 
Perhaps surprisingly, the larger the PSF (for the same galaxy 
size) the better the correction: presumably this is a small sam- 
pling effect. 

Finally, we introduced a lopsided PSF (by giving 1/4 of the 
photons an extra offset of half the FWHM) and repeated the 
analysis. As can be seen from Fig. [7] this combination of dipole 
and quadrupole PSF distortion can also be handled. 

In all cases, taking the expansion loN - 12 increases accu- 
racy, though not spectacularly. We conclude that the algorithm 
works: shapelets provide a promising technique for measuring 



galaxy ellipticities, and for correcting ellipticities for smearing 
by the PSF. 

4.4. Noise 

A final series of tests was used to check the behaviour of the 
algorithm on noisy images. We added Poisson noise to simula- 
tions such as those just described, and compared the result of 
many realizations. The PSF FWHM ranged from 4 to 12 pix- 
els, the galaxies' effective radii were set to 4 pixels. The same 
ranges of Sersic (0.5^-) and Moffat (2, 3, and 9) indices were 
used as above, and the input shear was 0. 1 . Different levels of 
noise were added, roughly to span the S/N range between 10 
and 100. 

The results are summarized in Fig [S] and [5] which show 
that the noise causes a scatter on the ellipticity estimates, but 
that this does not lead to a bias; and that the propagated error 
estimates on the ellipticities are a good measure for the rms 
scatter among the different realizations. 

5. The pipeline 

We have implemented the above ideas into a 'shear pipeline'. 
It starts from a fully reduced image, detects sources, deter- 
mines the PSF and its variation across the image, decomposes 
detected sources into shapelets, and obtains a shear estimate 
for each object. The pipeline consists of a number of stand- 
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Fig. 6. Residual shear after correction for an elliptical PSF (the same PSFs as fig.|U sheared by 10%). 
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Fig. 7. Residual shear after correction for a lopsided PSF (the same PSFs as fig. E] but a third of their flux is displaced by 
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Fig. 8. The result of Monte Carlo simulations, in which many noise realizations of Sersic profile galaxies were run through the 
ellipticity-fitting procedure described in this Paper. Shapelet order N = & was used throughout. Each plotted dot represents the 
average ellipticity of 2500 different noise realizations of the same galaxy image. The same data are plotted in both panels, but 
coded by different model parameters: the Sersic index on the left, and the PSF size on the right. The vertical axis shows the 
fractional scatter of the measured fluxes, cr(F)/F, of the sources, as determined by integrating their shapelet series. No trend of 
the mean shear with S/N is seen: noise leads to scatter but no bias. 



alone programmes that run in sequence. The modules operate 
on source catalogues and FITS images, and generate new cat- 
alogues and diagnostic plots. The pipeline can run fully auto- 
matically. 



5.1. Detecting sources 



The first step in the reduction process is the detection of 
sources. For this the SExtractor software (Bertin & Arnouts 
1996 1 is used. A few basic parameters are measured during ex- 
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Fig. 4. Illustration of the effect of the choice of scale radius 
on the ellipticity measurement. Each curve shows, for a differ- 
ent galaxy profile, how the derived e depends on the choice of 
scale radius (expressed as a multiple of the dispersion of the 
best-fitting round Gaussian). Each model galaxy had an effec- 
tive radius of 4 pixels and ellipticity 0.3, and was convolved 
with a PSF of Moffat index 2 and FWHM 8 pixels. The /3 for 
convolved source and PSF are both scaled by the same factor. 

traction: position, the flux and its error, the FWHM, 1 major and 

1 FWHM is determined by doubling the FLUX_RADIUS parameter 
of SExtractor, and not the FWHMJMAGE parameter which is prone 
to failure at low fluxes. 




0.001 L 1 1 — 1 — 1 1 — 1 

0.001 0.01 0.1 

Propagated pixel flux error 

Fig. 9. Comparison between the scatter in ellipticity measure- 
ments from sets of 2500 random noise realizations, and the er- 
ror predicted by propagating the pixel noise through the calcu- 
lations. 



minor axis length and position angle, and source quality flags. 
SExtractor is fast and effective, particularly for the relatively 
high S/N sources that are used here. 

5.2. PSF maps from shapelets 

The PSF is determined from the stars in the image itself. We 
assume that the pixel values are linearly related to intensity. 

First, the stellar locus is determined from the plot of mag- 
nitude vs. FWHM in the standard way (e.g., KSB). We fit a 
circular Gaussian to each star, and adopt the median dispersion 
from these fits, multiplied by 1.3 (see 34. 2> as the scale radius 

ySp S f. 

For all objects near the stellar locus a shapelet expansion is 
fitted using /3 ps f as scale radius. The SExtractor centroid is cho- 
sen as initial centroid for the expansion, but after the expansion 
is completed the centroid is adjusted (using the lst-order trans- 
lation operators of R03) until the coefficients of B 01 and B m 
are exactly zero. Finding the required offsets involves solving 
a simple linear equation. Each shapelet is finally normalized 
to unit integral, by analytically integrating the counts in each 
shapelet term and dividing by the total. 

Once this is completed, we have a shapelet description of 
(candidate) PSF objects scattered over the image. A map of the 
PSF variation across the image can now be made by interpolat- 
ing the shapelets. We have found that a straightforward polyno- 
mial fit, coefficient by coefficient, works well, though complex 
PSF variations may require more sophisticated schemes such 
as weighted nearest neighbour averages (Christen 2006), Pade 
interpolants (Hoekstra 2004}. or even physically-motivated 
model fits (Jarvis & Jain 2004 Jain, Jarvis & Bernstein 2005 1. 
During the fitting of the variation of each coefficient over the 
image, deviant points can be rejected, leaving a cleaned sample 
of PSF objects. 

The result of this step is a recipe for the shapelet coefficients 
of the PSF at any point in the image. 

5.3. Shapelet encoding 

Once the PSF is determined, all other detected sources are ex- 
pressed as shapelets as well. As for the PSF objects, a shapelet 
expansion centered on the SExtractor coordinates, is fitted di- 
rectly to the observed pixel values. The statistical errors on the 
pixel values are propagated through the least-squares fitting, 
leading to errors (and if desired, covariances) on the shapelet 
coefficients. In the case of well-resolved shapelets and uniform 
noise level across the source, the shapelet normalization is such 
that the rms error on each coefficient is the same, and the cor- 
relation between errors on different coefficients small. 

Each source is encoded into shapelets with scale radius de- 
rived as described in 34.21 All pixel values within a radius of 4/? 
from the SExtractor centroid (at least 10 pixels) are used in the 
fit. For efficiency reasons in the shear estimation step, the al- 
lowed /3 values are quantized: allowed values are /3 — 2"^/3 ps f, 

n = 0, 1 , 2, After fitting, the center of expansion for each 

object is shifted in the image plane by means of the shapelet 



10 



Konrad Kuijken: Shears from shapelets 



translation operators until the 01 and 10 coefficients are zero, 
as before. 

This procedure describes the source as seen in the image 
plane, i.e., after it has been convolved with the PSF and pix- 
elated. An alternative approach also make sense: to convolve 
all basis functions with the pixelated PSF, and fit the observed 
source image as a combination of those. This yields a shapelet 
description for the intrinsic, pre-seeing, object shapes (Massey 
& Refregier 2005 1. As long as the same procedure is followed 
for the sources and PSF, the end result should be the same: the 
deconvolved image will be free of pixelation and PSF. We pre- 
fer our approach because it leaves the covariance between the 
fitted shapelet coefficients small. 

All sources are encoded to the same shapelet order as the 
PSF (typically 8 or 12), in order to avoid signal-to-noise de- 
pendent smoothing effects. For faint sources, the higher-order 
coefficients will therefore be very noisy, but still unbiased. 

5.4. Shears from Shapelets 

With a description of the shape of each source and the corre- 
sponding PSF, the next task is to determine the intrinsic shape 
parameters that are needed for a weak lensing analysis. 

As explained in Section l3~3l we derive the ellipticity as the 
shear that needs to be applied to a suitable round source in 
order to fit the observed image. The fit is applied completely 

in shapelet space. We use a scale radius equal to (p 2 - /3 2 s{ ^j 
for the intrinsic, deconvolved, circular model galaxy, and nor- 
malize all sources to unit flux before fitting so that the c„ can 
be used as radial profile shape parameters. Only sources with 
f5 > jSp S f are used. The C„„,i convolution coefficients need to be 
evaluated only once per value of /3. 

5.5. Cleaning the catalogue 

The resulting catalogue of sources with shape estimators needs 
to be cleaned in order to remove sources which are affected by 
neighbours, edge effects, poor fits, etc. We apply various cuts: 

1. SExtractor Flags We first exclude all objects for which 
SExtractor raised a flag (due to neighbours, being close to 
the edge, saturation, etc.). 

2. Unresolved and faint objects Next size and S/N cuts are 
applied to the catalogue: typically all objects with best- 
fit Gaussian radius smaller than 1.1 times that of the 
PSF, and those with flux less than 10 times the flux error 
(as measured by SExtractor parameters FLUX_AUTO and 
FLUXERR_AUTO), are removed. 

3. Shape cuts For the next cut, for each source the fraction of 
power F „ at each order n in the shapelet expansion is cal- 
culated. An unusually high amount of power at high order, 
particularly for odd n, indicates a source whose shapelet 
expansion is affected by a neighbour. Thresholds are set for 
each F„ based on the properties of the ensemble popula- 
tion, above which sources are rejected. Typical values are 
F 3 > 0.05, F A > 0.2, F 5 >0.l,F 6 > 0.2. By construction 
F\ = 0, but to filter out peculiarly lopsided sources a max- 
imum can be imposed on the distance by which the center 



of the shapelet expansion had to be moved in order to set 
the 10 and 01 coefficients to zero. No cuts are applied to Fi 
since that would be similar to a cut on image ellipticity. 
4. Radial profile cuts As a by-product of the shear estimation, 
radial profile parameters cq . . . c^_2 are generated. These 
can be used to further excise peculiar objects. If the shapelet 
scales are chosen properly and the shear fit worked well, 
most of the flux of the pre-seeing, pre-shear model source 
should be contained in the C° term. Catastrophic failures 
of the fit, or problems with the setting of the shapelet 
scale, can be identified as peculiar values of co. Requiring 
|cq — 1| < 0.5 filters out such cases. 

5.6. Average shear 

Once individual estimates have been obtained for each source 
these need to be combined in some way to generate a shear 
estimate. 

Conceptually the simplest methods are (i) to average the e, 
and divide by (1 - (e 2 )) (eq. |9j, and (ii) to identify the mode 
of the ellipticity distribution (provided it is centrally peaked), 
which identifies the intrinsically round galaxies. 

A better technique is to form a weighted mean, where the 
weight is driven by the amount of information about the shear 
field each source provides. The scatter in the ellipticity mea- 
surements of sources is due to the intrinsic dispersion in shapes 
s e , and to measurement errors. The latter can be estimated by 
propagating the noise in each image through the fitting proce- 
dure, and the former can be estimated as the excess variance 
of e in the source population. We therefore adopt individual 
weights w - (s 2 + cr\ + cr\y l when forming the mean of all 
measured ellipticities. 

The same weighting is then applied to determine the value 
of (1 - (e 2 )) in eq.[9] So the shear estimate is determined as 



(s 2 e + o- 2 ]+ oir 1 



w(e\ + e 2 ) - w{cr\ + a- 2 ) 
(e > = — 



1 



we i 



V1Y'2 



1 - (e 2 ) w 



(18) 



(19) 



(20) 



The value of s 2 in eg. 1181c an be iteratively adjusted to equal (e 2 ) 
from eq.^5] though its precise value is of little consequence. 

Depending on the form of the intrinsic shape distribution 
of galaxies, different weightings are optimal: for example, for 
a very peaked distribution of ellipticities higher weight can be 
given to nearly round sources, whereas for a top-hat distribu- 
tion the sources with large ellipticity carry more information — 
for a discussion see BJ02. A problem with weight factors that 
depend on e is that the centroid needs to be found first as it is the 
intrinsic, pre-shear ellipticity that counts. The weight adopted 
above is optimal for a Gaussian distribution of intrinsic ellip- 
ticities. 

In the tests below we will compare the weighting scheme 
just described, and a simple unweighted median. To the extent 
that the median identifies the center of the distribution of el- 
lipticities of the source population, i.e., the intrinsically round 
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Fig. 10. Results from the STEP1 simulations, which are based 
on modeling of the optics of the CFHT. Each plotted point 
represents the average ellipticity of about 2200 sources in one 
STEP1 image. Shapelets to order A^ max = 8 were used. PSFs 
to 5 are, respectively, round, with coma, with astigmatism, with 
defocus, and with 3rd and 4th-order astigmatism, and results 
for images with applied shears of 0, 0.05 and 0.1 are shown 
as three clusters of points in each panel. The correction factor 
1 - (e 2 ) is about 0.93 in all cases. Results from cluster to clus- 
ter are not statistically independent, but within each cluster of 
points they are. 



sources, no (1 - (e 2 )) correction needs to be applied to the me- 
dian. 



Fig. 11. As Fig. [TO] but now the median ellipticity is used as 
shear estimator. The STEP1 galaxy ellipticity distribution is 
very peaked, which makes the median a very efficient estimator 
of the center of the distribution. 



6. Tests on STEP1 data 

A rather realistic test of the whole procedure is provided by 
the 'Shear Testing Programme' (STEP, Heymans et al. 2006 1. 
Phase 1 of STEP has produced a large set of realistic simulated 
images across which a constant shear (in gi) and PSF smearing 
has been applied. The PSFs mimic realistic optical aberrations. 
These images provide an important test, since unlike the tests 
presented in |4] the pre-seeing sources do not have a constant 
ellipticity with radius. Furthermore, the STEP images include 
photon noise, and the sources may overlap. The PSFs are also 
somewhat smaller (in pixels) than in the tests of 33.31 scale 
radii used range from 2.2 to 3.2. 
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Table 1. Summary of the results of the STEP1 simulations. For 
each PSF, the slope m\ and intercepts c, of the best-fit linear 
relation between input and recovered shear are shown. The c, 
have been multiplied by 100 for clarity. As no input g2 dis- 
tortion was applied in the STEP1 simulations, m-i cannot be 
measured. 



An analysis of the STEP1 data with an earlier version of 
this software was reported in Heymans et al. (2006 1. The main 
improvement in the implementation since then has been the use 
of polar shapelets in the shear determinations, which allows 
truncation effects to be curtailed properly, and the use of larger 
scale radii as discussed in 34.21 

Results of the use of the present pipeline on the STEP1 
data are shown in Fig. EH For each of three g\ shear values, 
and six different PSFs, 64 separate simulated images were run 
through the pipeline, each image yielding shear estimates based 
on about 2200 galaxies. Table^summarizes the results per PSF 
in terms of a multiplicative correction factor m, and an additive 
offset c. It can be seen that in general the method suffers from 
very little bias. For the non-elliptical PSFs (0,3,4,5), the recov- 
ery is perfect within the noise (m = 1, c = 0), which indicates 
that the correction for dilution by PSF smearing works. On the 
other hand, for the comatic (1) and elliptical PSF (2) there is 
a residual additive term, of nearly 0.003 and 0.005 in shear re- 
spectively. This result is consistent with that of the simulations 
in 34.31 In addition the elliptical PSF suffers from a multiplica- 
tive bias of nearly 4%. The origin of this discrepancy is not 
clear at the moment. 

The STEP1 analysis revealed that many methods show a 
small systematic trend between the error in the derived shear 
and the magnitudes and sizes of the objects (Heymans et 
al. |2)06). The technique presented in this Paper is no dif- 
ferent, as illustrated in Fig. ^] These trends are still a puz- 
zle, but it is clearly important to trace their origin and fur- 
ther improve the accuracy of the shear measurement.. Possible 
causes include ellipticity-dependent incompleteness in the cat- 
alogues, problems estimating the intrinsic ellipticity dispersion 
(e 2 ), magnitude- or size-dependent neighbour contamination, 
or residual systematic issues in the method itself. Further work 
addressing these issues is on-going. 

7. Comparison with other techniques 

The method described here has several advantages. It goes be- 
yond the KSB description of PSF anisotropy as a convolution 
with a very compact PSF, and is in principle applicable to all 
PSF shapes. The correction for the PSF is performed in a sin- 
gle step, which avoids the need to separate the effect of the 
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Fig. 12. The recovered shear for the STEP1 simulations for 
PSF0 and input shears of 0, 0.05 and 0.1, split up by source 
brightness (top) and size (bottom). In both cases there is a sys- 
tematic, so far unexplained trend. The upper panel in each plot 
shows the ellipticity dispersion correction factor derived for 
and applied to each bin. 

PSF into an anisotropic part that shifts the ellipticities, and a 
round part that causes a dilution of the source ellipticities. The 
forward-fitting approach of a PSF-convolved model for the in- 
trinsic galaxy shapes to the observations allows error propaga- 
tion. The fact that the ellipticities derived are 'geometric' in the 
sense of BJ02 (i.e., they represent the shear to apply to a round 
object in order to fit the source) means that there is no need 
to derive a higher-order 'shear polarizability', but instead the 
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response of an ensemble of sources to shear can be predicted 
simply from the dispersion in ellipticities. 

The use of shapelets in this method is not essential, but does 
help to speed up the calculations, and gives a natural framework 
for isolating the m — 2 components that carry the ellipticity in- 
formation about a galaxy. It also allows the source lists to be 
filtered efficiently based on objective shape criteria, and gives a 
robust way of interpolating the PSF shape across an image. In 
cases where the PSF cannot we well-described by a truncated 
shapelet expansion (for example, poorly-sampled space-based 
observations) it is possible to extend the technique, by perform- 
ing the PSF convolution of eq.^]in pixel space instead of in 
shapelet space (Kuijken 1999 1. 

As was shown by the STEP1 project (Heymans et al. 2006 1 
a variety of techniques can be used to derive shears from 
ground-based image data, with residual errors around the 1% 
level. While the focus of STEP1 was on different variations of 
the KSB method, we have shown here that the shapelets tech- 
nique can do as well. As KSB is expected to hit a fundamental 
level of systematic error (because of its good but imperfect de- 
scription of the PSF effects), which may well be inadequate for 
the next generation of weak shear surveys, it is worthwhile to 
look to higher-order methods. 

The approach we have taken here differs in several ways 
from that of Refregier & Bacon ( 2003 1 and Massey & Refregier 
( 2005 1, which is also shapelets-based. We perform the shapelet 
expansions to a fixed order, and prefer not to introduce signal- 
to-noise dependent thresholds that may lead to biases in the 
derived shears (S/N dependent truncation and averaging do not 
commute). Our non-iterative procedure is also much faster. Our 
shapelet expansions describe the observed, post-seeing images, 
which means that the coefficients are statistically almost in- 
dependent (whereas shapelet coefficients of the deconvolved 
shapes are necessarily correlated), ensuring that they are vir- 
tually unaffected by truncation of the series. We explicitly al- 
low for errors in the centroiding of the sources in our ellipticity 
estimates, and forward-fit the observed images to account for 
PSF effects rather than deconvolving them using a truncated 
PSF shapelet expansion. Finally, rather than modeling the re- 
sponse of a (model or empirical) population of images to a 
shear, we derive 'geometric' ellipticities (BJ02) for each in- 
dividual source, and average these only at the end. 

Compared to the Bernstein and Jarvis (BJ02) technique, the 
two main differences are the way ellipticity is measured, and 
the way PSF effects are handled. BJ02 derive geometric ellip- 
ticities by fitting an elliptical Gaussian form to the observed im- 
ages. This is similar to, but not quite the same as, a low-order 
shapelet fit as used in this paper. For sources with constant- 
ellipticity isophotes the BJ02 method is exact, whereas ours is 
accurate only to o(e 2 ) (see Fig.13. The benefit is a much faster 
numerical convergence. BJ02 correct the ellipticities for PSF 
convolution by means of a description of the PSF as a Gauss- 
Laguerre expansion (equivalent to polar shapelets). Optionally, 
the images can be convolved with a rounding kernel before the 
ellipticities and PSF are measured, in order to improve the ac- 
curacy of the PSF model. Thus the main difference with the 
approach here is that BJ02 first derive a post-seeing ellipticity, 
which is then corrected for the PSF; here we forward-fit the in- 



trinsic ellipticity in one step. BJ02 separate the effects of PSF 
anisotropy (ellipticity bias) and circularly-symmetric smearing 
(ellipticity dilution). They correct for the latter by assuming an 
intrinsic light profile of the source that is a perturbation about 
a Gaussian, expanded up to the kurtosis. The fact that here we 
model the full intrinsic radial profile of the source as higher- 
order shapelets should provide higher accuracy. 

The Kaiser (2000) technique is a more sophisticated kernel 
convolution technique, in which a convolution kernel is con- 
structed which turns an image into one for which the effect of 
pre-seeing shear on all sources is known exactly. This allows 
one to find the shear that makes the source ellipticity distri- 
bution isotropic — this is then the opposite of the amount by 
which the population was sheared on its way to the telescope. 
The method is in principle exact, and operates in pixel space. It 
appears to have received relatively little use thus far (Wilson et 
al. l20()lall200Tbl Dahle et al. EffiHl . 

8. Conclusions 

Shapelets provide a neat framework in which to describe the 
transformations that an astronomical source image undergoes 
until it is registered on a detector. Gravitational shear, convolu- 
tion with a PSF, and pixelation can all be modeled within the 
shapelets formalism. All these elements can be combined into 
an efficient algorithm for extracting image ellipticities that can 
be used for accurate gravitational lensing shear measurements. 

The implementation of these techniques into a working 
pipeline is presented in this paper. Tests show that the pipeline 
is able to recover input gravitational shears with very small cal- 
ibration error (of the order of a percent) and PSF residual (bet- 
ter than a factor of 30 in PSF ellipticity). 

It remains to be seen to what extent this approach can be 
applied successfully to diffraction-limited PSFs, which cannot 
be described easily with a shapelet expansion. A different set of 
basis functions for the expansion might be the answer. Further 
possible improvements are also under investigation. 
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Appendix A: Circular Shapelets 

The Cartesian shapelets S ab at order n — a + b can be written 
as inhomogeneous polynomials in x and y of combined order n 
times a circular Gaussian. The leading-order term for S ab (x,y) 
is (R03) 

2 n kx a y b e- r2/ ^ (A.l) 

where k is a constant that is independent of order, and r 2 = 
x 2 + y 2 . The circular shapelet at (even) order n is the unique 
linear combination of the S ah with a + b — n that depends only 
on r. to leading order in r 

C» =VZ^.("J.?\s l *-'(x,y) 



i/2) 

= 2 n kk'(x 2 + y 2)"/2 e -r 1 K2f} 1 )^ (A.2) 

where k' is chosen to normalize C" to unit integral over the xy 
plane. From the fact that rotation of Cartesian shapelets only 
mixes terms of the same order n — a+b it follows that the linear 
combination in eq. IA.2l is circularly symmetric at all orders of 



